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ABSTRACT 

An efficient numerical technique for solving the equations 
resulting from finite difference analyses of fields governed 
by Poisson's equation is presented. The method is direct 
(noniterative) and the computer work required varies with 
the square of the order of the coefficient matrix. The 
computational work required varies with the cube of this 
order for standard inversion techniques, e.g., Gaussian 
elimination, Jordan, Doolittle, etc. 

INTRODUCTION AND BACKGROUND 

General 

"Thermal spreading resistance" is defined as the conduc- 
tive thermal resistance between a source region and a sink re- 
gion in a solid where the geometry is such as lo preclude one 
dimensional heat flow. Knowledge of thermal spreading resist- 
ance is needed in two aerospace engineering areas. These are 
the thermal design of electronic components or equipments and 
in the prediction and control of thermal contact resistance. 

Importance to the Design of Electronic Components and Equip- 
ments- The thermal analysis of a power semiconductor or inte- 
grated circuit can be reduced to the problem of determining the 
appropriate spreading and bonding thermal resistances. As an 
example, the problem of calculating the junction- to- case ther- 
mal resistance of a semiconductor bonded to a substrate which 
is bonded in a metal case will be considered. Figure 1 illus- 
trates this problem. 

Heat is generated in a region of known size, the junction 
region of the semiconductor. The first, and most significant. 


*The work reported in this paper was done under National Aero- 
nautics and Space Administration contract NAS8-28616. 
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Fig. 1 Semiconductor in an Integrated Circuit 

spreading resistance of interest occurs between the junction and 
the opposite face of the silicon chip. The next thermal resist- 
ance of interest is that across the bond between chip and sub- 
strate. It is of significance that these thermal resistances are 
not independent. The thermal conductance of the bond proper 
can vary several thousandfold depending on the use of a metallic 
or nonmetallic bonding material. The resistance to heat flow be- 
tween the semiconductor chip bond region and the rear of the sub- 
strate represents a second spreading resistance. In a typical 
integrated circuit package, the entire bottom region of the sub- 
strate would not be available as a sink for a single semiconductor 
chip due to the presence of other heat dissipating chips. It is 
usually possible to estimate the size of the effective sink region 
on the rear of the substrate from considerations of symmetry or 
because it exceeds dimensions which appreciably affect the ther- 
mal spreading resistance. In those few cases where interactions 
must be considered, the key analytical tool is superposition. 
Green's function approach may be employed to advantage. For 
example, see reference 1. 

The importance of being able to predict thermal spreading 
resistances in single and multilayered material in the evaluation 
of the thermal design of semiconductor or integrated circuits has 
been shown. Spreading thermal resistances are important in 
other electrical devices such as phased array antenna elements, 
Peltier coolers, Seebeck generators, and most devices which 
utilize conductive heat transfer. 

Prediction and Control of Thermal Contact Resistance 

The resistance to heat flow between two mating (touching, 
as in a joint) pieces of metal is called thermal contact resistance. 
When the actual microscopic regions of contact between two mat- 
ing surfaces are examined, it is found that metal-to- metal con- 
tact occurs in small discrete regions where the asperities or 
microscopic protuberances make contact. References 2 and 3 de- 
scribe this model of contact in great detail. Figure 2 illustrates 
this contact model. 
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REOIONS OF ASPERITIC CONTACT 



Fig. 2. Microscopic View of the Joint of Contacting 
Pieces of Metal 


The heat flow to and from a region of asperitic contact is 
seen to be of the "spreading” type, hi fact, the effective thermal 
contact resistance of any contact may be considered as the sum 
of the parallel microscopic spreading resistances in the contacts 
themselves. References 2 arid 3 deal largely with isotropic con- 
tacts in which the thermal conductivity within the bodies of both 
contacts is uniform. 

Analysis has shown that the bulk of the spreading resist- 
ance occurs close to the region of actual asperitic contact and 
that the spreading resistance in any region varies inversely with 
the thermal ccxiductivity of the material. Figures 3 and 4 illus- 
trate the first of these points. Figure 3, drawn to scale, shows 
the equipotemial lines about a circular contact region eachdrawn 
to show one-tenth of the total spreading resistance between the 
circular contact region and a body of semi- infinite extent. It is 
seen that half of this total resistance occurs within one contact 
radius from the circular contact or source region and 80 percent 
occurs within three contact radii. Figure 4 illustrates these re- 
lationships. Figures 3 and 4 are taken from reference 4. 

The thermal conductivity of the contact material close to 
the contact region is of such importance that even a thin 45 Ang- 
strom thick layer of oxide on an aluminum contact can contribute 
measurably to the thermal contact resistance. This has been 
shown by actual measurement; see reference 4. 

Mikic and Carnasciali, reference 5, have utilized the above 
facts to enhance thermal contact conductance by plating materials 
of higher conductivity on the contacting faces of a metallic joint. 
They present an approximate analysis of spreading resistance 
from a circular contact into a contact region composed of two 
layers of materials with different conductivities. An exact bound- 
ary value solution of this basic problem has proven impossible as 
no mathematical function has been found which will satisfy all the 
boundary conditions between the plating and the body materials. 


♦ 
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Fig. 3. Temperature Profiles Described by Holm's Equation 
for Isothermal Circular Source on a Semi- Infinite Slab 

Professor C.J. Moore, Jr.,'*' in his discussion at the end 
of reference 5, felt this two layered spreading resistance problem 

'^Associate Professor of Mechanical and Aerospace Engineering, 
North Carolina State University, Raleigh, N.C. 
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DISTANCE FROM SURnCC- CONTACT RADII 

Fig. 4. Percent of Total Constriction Resistance for a Single 
Isothermal Circular Source on a Semi-Infinite Slab as a 
Function of Distance into Body of Contact 

could best be handled by a "well-conditioned finite difference 
computer code. " Mikic and Carnasciali then question the eco- 
nomic feasibility of such calculati(»is. 

Attempts by the author of this study to solve the two layer 
thermal spreading resistance problem using a finite differcuce 
approach utilizing Gauss-Seidel Iteration have shown the cost of 
digital computer calculation to be excessive for large nodal 
systems. 

METHOD 

General 

The governing differential equation for the thermal spread 
ing resistance problem is Poisson's equation. For those spread 
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ing resistance problems that are two-dimensional* or may be 
reduced to two-dimensional problems, the equation is: 


a^T _ QM, 


( 1 ) 


Consider a rectangular field subdivided into rectangular 
subregions as illustrated in Figure 5. The heat balance equation 
describing the heat flow among element m, n and its four princi- 
pal neighbors is; 


^^m,n ' ^m, n+l^^m, n ^ ^^m,n " ^m-l,n)V 

' ' m,n 

^^m, n ~ ^m,n-l^^m,n-l ^^m,n ” ^m+l,n^'^m+l,n 


= Q 


m,n 


( 2 ) 


where 


T is temperature 

Q'" is heat generated per unit volume 
X, y, z are spatial coordinates 

H, V are horizontal and vertical conductances, respectively 
m is row index 


n is column index 


„ is heat generated in node m.n 
in, n 


The convention for the horizontal and vertical conductances used 
is shown in Figure 6. 

Each of tue following observations will be helpful in under- 
standing the discussion which follows: 

(1) When any temperature Tm n known (e.g. , as a 

boundaiy condition), it will affect equation m,n by 

yielding a virtual heat generation term Qin, n^^kich 

is subtracted from the right-hand side of equation (2) 

where Q' „ is: 
m, n 


Q' = T (H -i V + i» 
^m,n m,n' m,n m, n 


“+■ V ) 

m,n-l m+l,n' 


( 3 ) 


*Thv method developed is applicable to three-dimensional prob- 
lems as discussed in this last part of this section. 
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Figure 5. Rectangular Field Divided into 25 Finite Elements 



Figure 6. Nomenclature for Nodal Interconductances 

(2) If the original physical field is divided into nodes of M 
rows and N columns, then: 

(a) There will be N*M linear equations. 

(b) There will be not more than N*M unknowns (fewer 
if some temperatures are initially prescribed). 
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(c) There can be as many different and distinct nodal 
conductances as there are interconnections be- 
tween adjacent nodes. 

Now, if the system of linear finite difference equations is 
written in matrix form (using the nodes of Figure 5) from left to 
right, top row to bottom row, as in reading English, a coefficient 
matrix results that has a pattern characteristic for field prob- 
lems described by Poisson's or LaPlace's equations. This pat- 
tern is illustrated in Figure 7. 

It was noted by Karlqvist (reference 6) that the matrix in 
Figure 7 may be partit ion^ as show n. It can be seen that each 
of the submatrices is «/NM x «/nM and the size of the coefficient 
matrix is NM x NM where the original finite element matrix was 
N X M in size. 

Derivation of an Efficient Technique for Exact Solution (tf this 
System of Equaticxis 

The submatrices shown in Figure 7 are defined in Figure 
8. Expanding the partitioned matrices (Figure 8) into a system 
oi equations, having normalized each equation with respect to 
the diagonal element: 

Tj - AjTg = Cj 

“® 2^1 ■'■^2 " ^ 2^3 ' ^2 

■® 3’^2 "^3 ■ ^ 3*^4 “ ^3 

• • • 

• • • 

• • • 


the general equation has the form: 


■®i'^i-l ■'^'^i" ^i'^i+1 ' ^i 

( 4 ) 

The first equatiem can be solved for T^: 


Ti = Cj + AjT2 

(5) 

and the i-th for T^: 


Tf = Ci + AjTj^j + BjTj_ j 

(e: 


The goal is to find a recursion relationship built upon sue 
cessive substitutions, which provides a solution for the i-th 
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Fig. 8. System of Submatr^ces in Matrix Notation 


unknown in terms of the (i -t-l)>th. That is; 
'"i = 


(V 


Examining equation (5) for above, it can be seen that: 


Ai = Aj 


and B'j = Cj 


The equation for T2 is; 

T, » Cj . A2T3 . BjT, (8) 

which, when written in terms of the equation for Tj, becomes: 


Tj ■ [1 - 82*1]' S’'.- * [' - 82''1]'‘ l"2 * 82®1] 

where I is the unit matrix. 

The general coefficients found in this manner become; 
1-1 


(9) 


and 


A5 =[i-b,a;.,J- a, 




Therefore: 


Ti = ♦ b; 


(a square matrix) 
(a column matrix) 

( 10 ) 


The temperature matrices (columns) are found st.::^ng at 
theVNM row (or 5th row of Figure 8). 

as A’ - 0 (11) 


T = B’ 

vmn vwir 




Defining O as the order of the coefficient matrix (e.g. , 
Square matrix of Figure 7), the system cT equations has been 
solved by operating on 3 •/&- 2 submatrices, each of which is 
the square root of the size of the original NM x NM coefficient 


452 



matrix. Three inversions of these submatrices are re- 
quired. The total number of multiplications - a good measure 
of the computer effort required during solution - is: 

No. of Multiplications = 30^ + - O + (12) 

This may be compared against other direct methods (see ref. 7): 


Method 

Number of Multiplications 
Required during Solution 

Gaussian Elimination 


Jordan 

‘o 

Doolittle 

+o^ -^o 

Cholesky 

+|o^ +^o 


Cornock's method (ref. 8), a triangulation type, alsomakes 
use of the characteristic pattern of submatrices which results 
during a finite difference solution for fields described by 
Poisson's equation. When the field properties are homogeneous 
and isotropic, Cornock’s method is very powerful since only one 
of the abc • submatrk js of order need be inverted. How- 
ever, for the general lution of the nonhomogeneous field the 
number n.alti»/,lcations required is: 

- 20- 5 (13) 

Cornock’s method does not lend itself to ready general prt^ram- 
ming for matrices of variable size as does the method described 
in this paper. 

That equation (12) 'S indicative of the computer c.lort re- 
quired for solution has been substantiated in practice. Figure G 
shows the variation is cost realized in the solution of very large 
matrices using the method developed in tbio paper. The mea- 
sured slope of the line in Figure 9 is very close to 2. 
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ORDER OF COEFFICIENT MATRIX (0> 

Figure 9. Cost of Computation versus Size of Coefficient Matrix 


Applicability of Technique to Three-Dimensional Problems 

The technique discussed above is suited to the solution of 
field problems having three or more dimensions. Figure 10 il- 
lustrates the characteristic pattern of the coefficient matrix for 
a three-dimensional finite element array. It is seen that the 
submatrices are 'v5~x “v^in size for an Ox O coefficient 
matrix. The same block tridiagonal pattern of submatrices is 
seen to occur as in the two-dimensional case so the derivation 
above for the technique of solution for two-dimeMional matrices 
is still applicable. Since the submatrices are x '^/5'rather 
than'^/5‘x '</5‘in size, the technique is even more powerful for 
three-dimensional problems. The number of multiplications re- 
quired for solution of the thr^e-dimensional problem is a func- 
tion of 0^/3 as opposed to Cr for the two-dimensional array 
where O is the order of the coefficient matrix. 
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